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1. INTRODUCTION 

This report deals with some calibration problems of monitoring a 
multiple array underwater tracking range. The arrays in the system are of the 
short baseline type; each contains four sonar transducers placed rigidly at the 
corners of a cube in a manner that describes a Cartesian coordinate system in 
three dimensions. Figure 1 contains a diagram showing the structure for 
these arrays and indicating the numerous signals they may receive. It is 
slightly deceptive in that real ray paths are not straight lines. 

An array receives a distinctive signal from a synchronously timed pinger 
attached to the target vehicle. The differentials of the sound wavefront's 
times of arrival at the four hydrophones allow the computation of the 
azimuth and elevation angles (spherical coordinate longitude and latitude) of 
the normal to the wavefront at the origin of the local coordinate system. 
Then, assuming direct path propagation, one can ray trace using Snell's law, 
[1], starting with the aforementioned elevation angle and utilizing a velocity- 
versus-depth profile for the speed of sound in the water. Finally, the time 
differential between the source pulse at the target vehicle and its arrival at the 
array is used to stop the ray-tracing algorithm and determine the location of 
the target relative to the array. The local track is the sequential set of these 
estimated positions. 

Each array in the system operates over a limited radius. As the target 
sojourns through the range, it is tracked by a number of these arrays. See 
Figure 2 for a plan view of the Nanoose range. (The zero level in the vertical 
is taken as mean sea level.) The overall path is constructed by transforming 
each piece of local track to the coordinates of the range based upon the 



assumed location and orientation of the various local coordinate systems. 
Discontinuities, or mismatches occur because the track produced by one array 
does not mesh well with that produced by a neighboring array in the overlap 
regions. Such mismatches can be seen in Appendix B. 

It appears that there are several sources of systematic error in the 
operation of this system. The question of individual array location and 
orientation has been treated previously [2]. The present report addresses the 
timing synchronization problem, from a statistical point of view. That is, the 
pulses received at the arrays are timed to the range clock with great precision. 
The timing information must be transmitted to the pinger prior to the release 
of the target vehicle. Although there is no engineering reason to suspect 
noticeable error in this transfer, some data exhibit behavior consonant with 
such an interpretation. Thus we build a mathematical model to account for 
such a source of systematic error, and develop statistical methodology for 
interpreting the data in the light of the model. Indeed, if the method 
provided a way to eliminate mismatches for an entire run for a single 
vehicle, there would be great temptation to use it as a smoothing filter. 

Generally there are several sequences of time points, called point count 
sets, for which two arrays simultaneously produce track. These occur for 
tracking in the overlap regions indicated in Figure 2. The paired tracking data 
of these point counts are called "crossover data." It is the crossover data 
generated as a result of the target vehicle's entire trip that provide the 
evidence suggesting timing synchronization problems and the data base for 
evaluating the use of such a model. Each version of track in a crossover data 
set is assumed to have been converted to the (common) range coordinate 
system. The components of this system are called "downrange, crossrange, 
vertical," or sometimes "centerline, crossline, vertical." 
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Figure 2. Nanoose Range 
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The model is developed in Section 2. It allows for both a timing offset 
and a drift. Estimation methods and statistical properties are developed. 
Section 3 contains applications to real data. It is shown that the offsets are 
significant and the drifts are negligible. The issue of the reality of the effects is 
also treated. I.e., is it the timing offset constant for all crossover regions used 
in a given run? To answer this question, some modeling of the components 
of variance is required and when done, a statistical analysis is performed. It 
appears that the true cause of these effects is something other than a timing 
synchronization bias. The analyses are supported by some characteristics of 
the noise process, and these are presented in Section 4. Conclusions are 
summarized in Section 5. A number of appendices are included that contain 
supporting data and information, including source codes for computations. 

2. MODEL DEVELOPMENT AND STATISTICAL PROPERTIES 

Figure 3 contains a mockup illustrating conditions that support the 
consideration of a timing synchronization correction. It shows a plan view 
which includes the radial lines from the arrays to the estimated track points 
in adjacent overlap regions. Between these regions, much data is supplied 
only by array A 2 . The analyst does not see the radial lines on his screen, nor 
does he see the black dots. He sees only the X's and O's (no distinction 
between the two) and no mismatch is apparent. (Mismatches would be 
apparent however for track pointing in a different direction.) But when one 
pairs up the radial lines by common time points, then one sees that the two 
versions of track lack coherence and can be improved by stretching the 
estimated points to the black dot positions. This can be achieved by a single 
constant adjustment to the transit time values in the ray tracing algorithm. 
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Figure 3 : Effect of Timing OffseL Errors (Synthesized) 
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Track according lo Arrays Aj or A3 as indicaled 



The figure also helps one to imagine the distinction between the effects of 
an array location error model and a timing adjustment model. If the former 
were applied to the situation in Figure 3, at least two arrays would have to be 
moved. But anv such decisions cannot be made in isolation. The ultimate 
goal is the simultaneous improvement of coherence in all overlap areas and, 
of course, discontinuities must not be created at other places. 

A crossover data set is a set of matched pairs of three vectors 



fXi(t)l 




Y,(t)' 


X(t) = X2(t) 


^ Y(t) = ' 


Y2(t) ^ 


[X3(t)J 




lY3(t)J 



representing two versions of the same vehicle track for a common set of time 
values (point counts), T in number. The array that produced X(t) is located at 
a and the one that produced Y(t) is located at p. Thus these two versions of 
track can be represented in their local coordinate systems as 

^(t) = X(t)-a r\it) = Y(t)-P (2) 



Timing synchronization error corrections may be viewed as either stretching 
or contracting the vectors ^(t) and ri(t) by the same (time) amount. The effect 
of such corrections is not constant, but depends upon (i) the speed of sound in 
water at the depths X 3 (t) and Ysft); (ii) the elevation angles of the ray traces at 
these times and depths. 

The magnitudes of the time adjustments are small and the effects can be 
described using first order terms 



C(t,g(t)) = c(t) + g(t)a(t) 
Tl(t/g(t)) = ri(t) + g(t)b(t) 



(3) 



where g(t) is a scalar function of time representing the total adjustment for 
timing offset (5) and drift (m) 



g(t) = 5 + m(t-to) (4) 

for a conveniently chosen central time, to; and a(t), b(t) are the scaled 



directions of stretch 



a(t) = v(t) 



cos(6(t))cos«)(t)) 

cos(0(t))sin((t)(t)) 

^sin(0(t)) 



(5) 



i.e., v(t) is the speed of sound at depth ^3(0, 0(t) is the elevation angle of the 
ray trace from a at vs(t), and 0(t) is the azimuth of w(t) from a. The vector b(t) 
is defined similarly for the track ri(t) and its origin p. 

It is easily seen that, for given 6, m and to, the corrected versions of track 
in the (common) range coordinate system are 

X(t,g(t)) = X(t) + g(t)a(t) 

Y(t,g(t)) = V(t) + g(t)b(t) 

The statistical estimation problem is to chose 5 and m so that the two 
versions of track agree as well as possible. The least squares approach is 
adopted. Using the notation 

I ^ 

A''e|:lV(()[f = -In''(0VV(0 C) 

' ' /=! 

we set up the objective function 
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( 8 ) 



Q = Ave|x(/,^(/))-y(/,g(/))f 
= Ave||X(/)-y(/)f + Ave[g(/))^K/)-KOf 

+2Aye x(0[«(0-K0] [>^(0-'^(0l 

and find the values of 5 and m that minimize Q. When the two partial 
derivatives are set equal to zero it is seen that there is convenience in 
choosing to so that 

Aye(t-to)||fl(/)-h(Of = 0 (9) 

That is 

to = Aye t\\a{t)-b{t)f / Ave ||a(t)-h(/)f 

This done, the normal equations may be expressed as 

5Aye||i7(f)-h(/)f = - Aye(a(/)- h(/)) [X(0-V(/)] 

( 10 ) 

m Ave(t - /o)‘|KO - KOf = - Aye(/ - fo)[^(0 " MO] ^'(0] 

and one can readily solve explicitly for the minimizing values S and m. 
Positive values of S are associated with stretching, negative with contracting. 

If we let Q be the minimizing value of Q, it is useful to establish the 
decomposition 

Q = Q + (^ - 5) Aye||fl(/) - b{t)f + (w - Aye(/ - /q )]«(/) - ^(Of 
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Before justifying (11) it is convenient to shorten the writing: let U’(/,(5(f)) 
= X(/,^(/))- y(/,g(/)); c{t) = a{t) - b{t); g{t) = 5 + m[t - {q). Then the statement 



sa)'s 

Ave||w(/,g(/))||‘ 

= Ave]|lV(/,g(/))|j‘ + (<5 - <5)‘ Ave||r(f)f + {m - mf Ave{t - /o)|K0||‘ (^2) 

and the result is justified by use of the three orthogonality relationships 



Ave 

t 



lV(f,0) •+• 5c{t)-r m[i - /o)c(/) c(/) - 0 



Ave 

t 



lV(/,0)-r &-(/)+ - to)c{t)\ [t - to)c{t) = 0 

Ave(/-foy'(0^(0 = 0 



(13) 



which in turn are established using the normal equations (10), and (9). 

The significance of the offset and drift parameters can be judged if we 
develop the means and variances of 5, in. Letting E denote the mathematical 
expectation operator, we begin with the assumption 

E[W(t,g(t)] = 0 (14) 



which embraces the idea that the two corrected versions of track produce 
common track without systematic error. 

Further, let 



-1 
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rn 



1-1 






05 ) 



and use (13) and (14) to show that the estimators 5 and in are unbiased. 
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e(s) = -/C5X^'(OC[W(f,0)] = KsJ^c'itp + m(f - fo)]c(f) = <5 

t t 

E(A) = -K,„X((-<oy'(')£[‘V(f,0)] = K„£(/-(„)c'(r)[5 + m(l -/„)):(() = >n 

t t 



To develop variances, we assume that the positive lag covariances are 
zero; support for this appears in the section on noise characteristics. Let M 
represent the (zero lag) covariance matrix. (See Appendix C for estimates). 



M = £[w(Lg(f))lV'(/,^(0] 


(16) 


Then one easily represents 




var(^) = xj£c'(f)Mc(t) 
t 


07) 


var(m) = c'{t)Mc'{t) 

t 


(18) 


cov( S,m) = KsK^"£{t 


(19) 



( 



If M is proportional to the identity matrix then this last term is zero; the other 
two terms simplify immensely; and a standard regression development can 
be used. But the study of noise characteristics does not support this. On the 
other hand the vectors (c(t)} do not change much with t and this has the 
tendency to render (19) to be small. The reason for this stability is that 
crossover data occurs only at the greater distances from the sensing arrays; the 
azimuth and elevation angles and the layer sound speeds do not change 
much. 
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It appears that the matrix M changes with the crossover data set. 
Methodology for estimating it appears in Section 4. Estimates of the M 
matrices appear in Appendix C. 

3. DATA ANALYSIS 

The data consists of S, m , and supporting values for 69 segments of 
crossover track collected over four separate days with two (temporally serial, 
not concurrent) target vehicles (runs) per day. The information is 
summarized in Table 1. Missing variance estimates indicate either a data 
shortage or outlier problems. In a few cases the tracks were curved, and the 
straight line model is not adequate. The units are milliseconds for 5 , and 
milliseconds per point count for rh. One millisecond translates to about 4.5 to 
5 feet of distance. The ratios of means to standard deviations are used to 
judge whether the effects are significantly different from zero. Virtually all of 
the offsets and some of the drifts are significant although the latter are not 
strongly so. The "r” column contains the correlations between 5 and m . 
They are insignificant. A further search for large scale drift was made by 
plotting 5 against the crossover central time to for each run. They appear in 
Figure 4. If a smooth signal were discernible then we would have a way to 
connect the 6 values that appear in each column of Table 2. But they are 
chaotic and provide no incentive to continue any concern about drift. It is 
concluded that the offsets are significant and the drifts are not. The latter will 
be dropped from further consideration. Some graphical examples of the effect 
our timing corrections have been selected and appear in Appendix B. 
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Table 1: Offsets, Drifts, Signal to Noise Ratios 



Dale 


Vehicle 


6 






m 




m/Orr, 


r 


^0 


T 


3/23/89 


MK 30 


1 37 


0461 


29 82 


00168 


.001627 


1 03 


- 00 


405 1 3 


4 2 






1 74 


.0431 


40.28 


- 00831 


00221 7 


3 75 


.00 


422 72 


3 3 






1 69 


0925 


18 29 


- 01 559 


.007024 


2.22 


.06 


4171 09 


1 9 






78 


0464 


16 87 


00245 


001 588 


1 54 


02 


439 71 


45 






- 1 5 


.0504 


2 97 


.01304 


001438 


9 07 


02 


201 2 80 


57 






10 


.0802 


1 23 


00150 


.005509 


27 


01 


3080 16 


22 






1 24 


0480 


25 96 


- 001 1 9 


001645 


.72 


01 


2531 .21 


49 


3/23/8 9 


MK 30 


3 54 


.0479 


73 85 


- 00395 


0021 19 


1 .86 


06 


422 64 


3 7 






2 52 


0584 


43 17 


- 00224 


0021 36 


1 05 


1 1 


1 605.22 


4 7 


5/10/89 


MK 27 


3 04 


0283 


1 07 40 


00993 


001 1 42 


8 70 


♦ 02 


31 07 58 


4 3 






3 44 


0192 


1 78 88 


- 00250 


000866 


2 89 


- 06 


31 20 22 


38 






3 24 


.0401 


80 68 


• 00007 


.002747 


02 


.05 


6936 70 


25 






2 43 


0270 


90 18 


00187 


001375 


1 36 


00 


31 1 6 86 


34 






1 86 


0266 


69 80 


.00029 


000825 


.35 


-.00 


3781 .03 


56 






1 24 


0228 


54 56 


00408 


000630 


6 48 


- 05 


5471 41 


59 






2 06 


0197 


104 59 


.00437 


.000662 


6 60 


.02 


4720 56 


50 






1 99 


0104 


190 86 


- 00067 


000364 


1 84 


02 


8060 02 


50 






2 27 


0180 


126 16 


-.00196 


000770 


2 54 


- 03 


5992 35 


40 






1 89 


0868 


21 72 


- 00001 


000567 


.01 


- 01 


961 5 69 


50 


5/ 1 0/89 


MK 27 


3 49 


01 94 


180 1 1 


01062 


000746 


14 23 


- 03 


1733 39 


45 






3 66 


0270 


1 35 56 


- 00044 


001 268 


35 


00 


1751 02 


37 






4 90 


01 89 


259 57 


.00443 


000689 


6 43 


01 


10196 44 


4 5 






2 93 


0258 


1 13 27 


- 00545 


001358 


4 02 


00 


1 747 85 


33 






2 09 


0406 


51 50 


00038 


001887 


20 


• 01 


2366 96 


3 7 






2 85 


0243 


117.13 


- 00094 


000845 


111 


07 


2980 01 


50 






3 87 


.... 


..... 


- 00050 




- . - • 




1 2354 91 


50 






2 37 


.... 


..... 


- 01074 





- . - - 




5316 80 


1 5 






3 38 


0131 


257 71 


00257 


000458 


5 62 


00 


12215 15 


50 






1 36 


0209 


64 92 


.00276 


.000662 


4 1 7 


- 06 


3585 00 


55 






2 1 8 


0181 


1 20 55 


- 001 95 


000628 


3 1 1 


07 


4478 00 


50 






2 18 


0149 


1 46 48 


00192 


000479 


4.02 


-.03 


5873 92 


54 






4 91 


0096 


510 80 


00022 


000337 


64 


02 


11235 86 


50 






3 1 5 


0281 


11197 


- 00202 


001001 


2 02 


- .05 


6480 62 


46 


6/6/89 


MK 27 


2 68 


0336 


79 78 


01680 


003073 


5 47 


.02 


1952 64 


1 9 






4 31 


01 1 6 


371 .65 


- 001 89 


000377 


5 02 


08 


10562 53 


50 






4 95 


0269 


18419 


001 44 


.001 1 36 


1.27 


- .04 


1 441 1 90 


4 1 






2 72 


01 25 


21 7 40 


00343 


000434 


7 9 1 


02 


2124 12 


50 






4 03 


0136 


297 47 


00242 


000464 


5 22 


-.02 


1 2709 75 


50 






1 1 4 


0376 


30 32 


- 00197 


001 555 


1 27 


- . 01 


4979 22 


40 






1 91 


0251 


76 25 


.00401 


001141 


3 52 


07 


61 97 42 


3 5 






2 69 


01 51 


1 78 04 


- 001 43 


000486 


2 93 


- 05 


1 2044 01 


5 1 


6 '6 8 9 


MK 27 


2 1 1 


0232 


91 30 


00091 


000783 


1 17 


03 


357 26 


50 






3 09 


0434 


7115 


00596 


003137 


1 90 


-.02 


421 96 


24 






2 90 


0195 


1 48 86 


- 00231 


000677 


3 42 


- 03 


41 35 53 


50 






2 48 


0233 


106 78 


-.0021 1 


000692 


3.06 


04 


487 61 


50 






1 99 


.0156 


1 27 36 


- 00192 


.000541 


3 55 


- 02 


1781 82 


50 






1 .04 


0202 


51 69 


- 00075 


000681 


1 10 


.04 


2505 37 


50 






1.32 


0181 


72 74 


00003 


000631 


05 


01 


1991 06 


50 


7 ' 2 1 / 8 8 


MK 27 


2 07 


021 7 


95 45 


00060 


000725 


83 


- 02 


352 1 5 


52 






2 75 


0300 


91 58 


00020 


001845 


1 1 


- 00 


377 91 


28 






2 24 


01 1 5 


195 77 


- 001 60 


000398 


4 01 


- 00 


4440 01 


50 






1 90 


0147 


1 29 27 


00244 


000510 


4 79 


.03 


469 17 


50 






1 47 


.... 


. . . . 


- 00746 




- . - - 




1062 61 


1 5 






1 16 


0128 


90 47 


- 00287 


001310 


2 1 9 


01 


1 772 98 


1 7 






29 


01 44 


1 9 86 


00232 


000501 


4 64 


-.01 


2068 58 


50 






2 45 


01 1 2 


217 86 


00069 


000392 


1 76 


00 


3137 66 


50 


7 '21 /8 8 


MK 30 


1 37 


0245 


55 85 


01 1 49 


000853 


1 3 46 


02 


1108 15 


50 






2 1 5 


0306 


69 77 


.00069 


002324 


.30 


- 02 


1 1 52 98 


23 






96 


0192 


50 1 8 


00019 


000662 


29 


- 03 


6658 19 


50 






1 10 


01 71 


63 97 


00206 


000595 


3 46 


.02 


217 45 


50 






69 


0189 


36 48 


- 00233 


000656 


3 54 


-.02 


4360 14 


50 






02 


0099 


2 07 


.00190 


000346 


5 49 


- 06 


1 1498 22 


50 






- 80 


0209 


38 31 


00147 


000730 


2 01 


■ 00 


4667 48 


50 






- 42 


0337 


12 35 


00355 


001 1 73 


3 02 


- 03 


10430 20 


50 






27 


0176 


1 5 46 


00035 


00061 1 


57 


- 07 


3793 99 


50 






1 7 


0251 


6 69 


00024 


000952 


26 


03 


5141 88 


27 






54 


0496 


1 0 88 


001 71 


003582 


48 


03 


1 2031 94 


24 






61 


0259 


23 51 


00509 


000897 


5 67 


02 


301 1 04 


50 






47 


02^ 


1 5 94 


. 23 


001 380 


2 34 


04 


5697 10 


3 7 
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Table 2: Estimates of the Timing Offset 



4.91 (i) 
4.89 (k) 
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Figure 4. Offset vs. Central Point Count for the eight runs 
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Since the 5 values are significant, it is important to discover what they 
really represent. They have been labeled timing synchronization offset errors 
and, if that labeling is a valid physical description, then each run (vehicle) on 
each day should have its owti unique value and this common value could be 
used to correct all track, not just crossover track. What follows is an analysis 
of this point. 

Table 2 contains a graphical positioning of the sixty-nine 6 values 
according to the eight runs. There were eleven array pairs involved in this 
process and they are marked with the letters a, b, ..., k. (The correspondence 
with the overlap regions marked in Figure 2 is given in Table 3.) Study of 
Table 2 shows considerable "within run" variability, and one is inclined to 
doubt the reality of the constant offset interpretation. 



TABLE 3. IDENTIFICATION OF ARRAY PAIRS TO RECONCILE TABLE 2 

WITH FIGURE 2 



a 


0,11) 


e 


(3,4) 


i (5,56) 


b 


0,2) 


f 


(43) 


j (5,14) 


c 


(2,11) 


s 


(5,6) 


k (15,16) 


d 


(2,3) 


h 


(435) 





In passing we also note that the 5 values for each array pair appear to 
have some temporal coherence, and hence the assignable causes may be 
related to this, but that is an issue for another time. 



We proceed to model the components of random error variance and 
develop a test statistic for examining whether or not all 5 for the same run 
can be vievsed as constant. The {(5, ) are modeled as being affected by the day 
(water depth velocity profiles change from day to day), the run (the second 
run is a different vehicle operating later in the day), and the array pairs 
generating the crossover data. The array pairs are assumed to produce fixed 
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effects, with values a, for j=l, 11 ; and the day and run effects are assumed 

2 2 

to be random with zero means and variances O] and 02 respectively. There is 

2 

also an experimental (residual) error having zero mean and variance a . 

The result is the mixed model 

5 = Xa + Zp + e (20) 

where <5 is the 69 component column vector of 5 estimates; X is a 69 x 11 
matrix of zeros and ones (incidence matrix) relating the S component to the 
array pair; a is the 11 component vector of array pair fixed effects; Z is a 69 
rowed partitioned incidence matrix 

Z = (Zi, Z 2 ) (21) 

with Zi having 4 columns relating the 5 component to its day and Z 2 having 

If 

8 columns relating the 5 component to its run. The vector P' = (Pj, P,) 
represents the random variables 

P, = (Pn, ..., P 14 ) independent N( 0 ,O]) 

2 

P 2 = (P 21 / •••/ Pis) independent ^( 0 , 00 ) 

and the residuals are 

e' = (ei, ..., e 69 ) independent N(0 ,g-) 

Further the vectors Pi, P 2 , e are assumed to be mutually independent. 

We will be applying the model (20) separately to each of the eight runs for 

purposes of estimating the fixed effects (a) for each of those runs. Such 

2 2 

estimates will require values for the variance components (a^, ai, ap. We 
prefer to estimate them but once, by pooling all of the data. This can be 
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accomplished separately from the estimation of fixed effects by using the 
) restricted maximum likelihood method: the unequal block size version has 
I been treated by Patterson and Thompson, Ref. [4] and is used here. (See 
Appendix A for a description of the algorithm.) The results are 



and the units are seconds squared. It is noted that the run-to-run variance is 
double that of the day-to-day variance, which in turn is larger than the error 
variance by a ratio of about seven to five. But the degrees of freedom for the 
first two variances are small and the estimates mav not be verv reliable. 

^ j 

Using the model (20) one sees that the covariance matrix of the 
observables is 



where U is the identity matrix of order n, T is diagonal, the first four values 
2 ^2 

being and the next eight o?; and Z is partitioned (Zi, Z 2 ) as before. For the 

full data set n = 69, but recall that our goal is to check whether the fixed effects 
(offset estimates) can be viewed as constant within each run. The plan is to 
estimate the fixed effects for each run and test whether they can be viewed as 
constant (over the array pairs). 

To do this we proceed as follows. The 5 values for the run can be 
identified by the ones in the k'*^ column of Z 2 , k = 1, ..., 8. The array pairs 
involved (not necessarily all eleven) are identified by ones in the X matrix 
restricted to the k^*^ run. The Z matrix is restricted in a like fashion and the 
covariance matrix (of order n, the number of crossover data sets in the k''^ 
run) has the same structure as (23) with the reinterpretation of inputs. Of 
course the estimates (22) must also be inpiut. 



Oi = 0.3571 02 = 0.7207 o = 0.2430 




(22) 



H = o^in + zrz' 



(23) 
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Let K be the number of array pairs involved in the run. Use the 
Aitken estimator (3] 

(24) 



and its covariance matrix 



var 



(“»r) = 



x'lHk'Xt 



-1 



= V 



-1 



(25) 



and the subscript k modifies the previous definitions of symbols so that only 
the k*^ run is involved. If the null hypothesis is true (i.e., the all fall on 
the main diagonal of the K dimensional space), then the maximum 
likelihood estimate of the common value is 

i=l;=l : / 

It also follows from the normality assumptions that the quadratic 

-a) V[ccj^ -a) - Chi Square (K-1) (27) 

and this statistic is the weighted distance of the components of 6. from the 
main diagonal (i.e., the constant fixed effect that is the same for all array 
pairs). Thus the null hypothesis should be rejected when this distance is too 
large. 

The numerical results for the eight runs are contained in Table 4. The 
column marked p* contains the p values (empirical significance level) 
corresponding to the test statistics listed under the "distance" column. Seven 
of the eight values indicate rather rare events and the eighth, p’^ = 0.14, is 
associated with one degree of freedom. Low degrees of freedom tend to 
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dampen the chances of finding significant results. The final column is the 
ratio of the distance (27) to its standard deviation under the null hypothesis. 

We assert that the array pair effects are not constant and that there are 
other sources of systematic error that dominate this process. 



TABLE 4. TESTING THE EIGHT RUNS FOR CONSTANT FIXED EFFECTS 



a 


K-1 


Distance 




distance/SD 


0.61 


6 


30.17 


3.6x10-5 


8.71 


3.03 


1 


2.12 


0.14 


1.51 


2.35 


8 


17.71 


0.06 


4.43 


3.09 


9 


51.67 


5.2xl(H 


12.18 


3.05 


4 


41.48 


2.1x10-8 


14.67 


2.13 


6 


14.16 


0.03 


4.09 


1.79 


7 


18.13 


0.01 


4.85 


0.46 


8 


30.44 


1.8x10-^ 


7.61 



4. CHARACTERISTICS OF NOISE 

In order to study the stochastic nature of the tracking errors we fit straight 
lines to the track. (This was done even for track that did not appear linear via 
a visual scan of graphical output. The exceptional cases are marked in 
Appendix C.) It is assumed that the target vehicles speed is constant. The 
result is a set of deterministically spaced points that fall on a straight line in 
three space. From these we can produce residuals and, assuming local 
stationaritv, studv the noise structure. 
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The basic (unadjusted for constant speed) straight line is found by 
principal components. Let X(t) be used for X(t,0), eq. ( 6 ) and 

X=lxx(() ( 28 ) 

' t 

VVe seek a projection (direction) p = (pn pz- Pa)' 

Z(<) = i.P.XM = P'Ml) (29) 

1 

that will maximize the variance of (Z(t)} 

<^i = 7l[z(0-zf = P'CxP (30) 

‘ t 

where 

Cx=7Xb’(0-X)(X(/)-X) (31) 

‘ t 

is the covariance matrix of the track data, X(t). 

Since the vector p is a set of direction numbers, we adopt the usual 
constraint, = h (This will keep finite.) It is easily shown that the 

three solutions to the eigen problem 

CxP = 

provide us with an orthonormal basis 

P = {Pl.P2.h] (32) 

2 

where P, is the eigen vector corresponding to />, and X] > X 2 > X 3 . Since 02 = 

2 

p'Cxp = Xp'p = X, and our goal is to maximize 02 . vce choose the first eigen 
vector for use. 
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The set of values 



Z,(0=P,'(X{1)-X) (33) 

represents the succession of values after the vectors |x(/)-X| are projected 

onto the line of the first principal component. These projections are 
orthogonal projections and may not well represent the relative positions of a 
target vehicle heading in the direction Pi at constant speed. Accordingly it is 
appropriate to make an adjustment in the {Zi(t)) to account for constant 
speed. The collection of point counts, 

t-i < t2 <■..< tj, (34) 

all differ by integral multiples of some base value, A, representing distance 
traversed per point count. Let us perform a simple linear regression of the 
Zi(t) on the times (34). The fitted values are 

Z{t) = a + b(t-t) (35) 

where = — X,Zi(f) and - r)Z](f) - f Xow the successive 

values of Z{t) differ by integral multiples of a converted base value, bA, and 
represent the constant speed progression of the vehicle in the direction P]. 

It remains to represent the distances Z(f) in the original coordinate 
system. Let 

R(;) = (z(/),o,o)' (36) 

be the estimated position of the vehicle in the basis P, and then 

X^{t) = X^PR{t) (37) 
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will be the straight line in the original coordinates. Similarly, when (Y(t)} is 
put into this algorithm, we produce Ysi(t). 

Now we are positioned to estimate the covariance matrices of the noise 
processes 



Computational work shows considerable variability in these matrices as one 
changes the day, the run, and the array pair. Also the cross correlations are 
mostly different from zero. The matrix M, eq. (16), is estimated by 



and used in (17), (18), and (19). The quantities (38) and (39) appear in 
Appendix C, and illustrate their variable nature. 

Another immediate use of the noise processes is to look at the 
autocorrelations. In addition to the three individual components we define a 
"noise displacement" process 









(38) 






M - + Dy - - Dy^ 



(39) 




(40) 



1 



(Of course the same is done for (Y(t)}.) 
The autocovariances 
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( 41 ) 



‘ f=l 

for {d(t)) and for the individual components {X,(t)-Xsi.(t)) are computed and 

normalized. Plots of R(h)/R(0) vs h appear in Appendix D. It is typical that 
they become and remain in a general level of "static" for h > 1. 

5. SUMMARY 

The paper contains methodology for estimating the presence of timing 
synchronization error and judging its significance in the framework of our 
short baseline array underwater test range. The methodology was applied to 
real data. Some plots illustrating the effects appear in Appendix B. One must 
view them with care because the coordinate scales are so different. It is rather 
typical that some unresolved systematic error is exposed in the side view. 

Analysis of the results does not support the timing s\mchronization error 
model as accounting for the discrepancies. Other sources of systematic e Tor 
must be unmasked first. Possibilities include the orientation of the arrays [2], 
biases in the raytracing inputs [3], and perhaps a temporal or spatial gradient 
in the depth velocity profile. 

Appendices C and D contain information about the second order 
properties of the three dimensional noise process and about time lag 
correlations. 
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APPENDIX A. PATTERSON AND THOMPSON ALGORITHM 



TT»e purpose of this appendix is to document the computational method 
for estimating variance components that is but implicitly described in the 
Patterson and Thompson paper [4]. We use the notation of that paper. 
Specifically 

y = Xa + e (A.l) 

where the observeable y is an n component column vector, X is an n by t 
matrix of rank t determined by the allocation of treatments to units, a is a t 
vector of fixed effects, and g is a mean zero n vector of normal random 
variables with covariance matrix 

V = o^H H = ZrZ' + 1. (A.2) 

The matrix I is the identity of order n, Z is the n by b design matrix for c block 
factors, and T is a diagonal matrix containing the variance components 
relative to the basic error variance a-. I.e., 

C 

H = 1 + X^pZp'Yp (A. 3) 

P=i 

and each block design matrix, Zp, is n by bp; each plot having exactly one level 
in each Zp, p = 1, c, the variance of the p'*^ block is 

Op = Vpo2 for p = 1, 2, ..., c (A.4) 

and Z in (A.2) has the partitioned form 

Z = (Zi, Z 2 , ..., Zc). (A. 5) 
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The diagonal elements of T are bi consecutive yi's, b 2 consecutive Y 2 's, and 

C 

be consecutive Yc'S/ with = b. 

P=i 

Out goal is the estimation of Y 1 / Y2/ •••/ Yo by the restricted maximum 
likelihood method for the unbalanced case. The steps follow. Let 

S = I-X(X'X)"’X' (A.6) 

W = Z'SZ+r"’ (A.7) 



The algorithm is an iterative one, and initial values for the Yi/ •••, Yc^re 
needed to develop W in (A.7). The b by b matrix W and its inverse can be 
partitioned according to bi, b 2 , ..., be. So doing allows both 





(i = VV ’z'Sy 


(A.8) 




u = 


(A.9) 


to be viewed as partitioned. 


I.e., 






P'= P,,P2 P, 

V ^ ^ 

D = {L'ij) i, j = l,...,c 


(A.IO) 


and L'ij is a bibyb. matrix. 






Now we are positioned 


to define 






R = y'Sy-y'SZp 
Pp “ Pp t^p/fp 


(A.ll) 



and Ep = trace (Upp) for p = L ..., c. Then define the modified information 
matrix (f,j}, of order c+1, 
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(A.12) 



f,j = trace (U,,Uji) for i, j = 1, c 
fp c+i ~ ^c+i/p ~ Ep for p — 1/ c 

^c+l,c+l ~ 



and let (f*)} be its inverse. 

Finally, the estimator update equations are, using k = c+1 




< for; = l,...,c 



(A.13) 



and the new {y,} can be placed in (A. 7) to start the next iteration. No initial 
value for is required, but it must be updated at each iteration. The 
algorithm should be stopped when (A.13) is stable. The variance components 
are computed from (A. 4) 

The application in the present report has c=2, with bi = 4 and b 2 = 8. It is 
helpful to record the partitioning and inversion of W, (A. 7). 



z,sz, z,sz, \h/r, 0 
Izjsz, zjsz; [o /j/ri 



(A,14) 




s}Tnmetric. 




(A.15) 



can be computed from 
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(A.16) 



W“^ = |vV22 — ^^12 1 

w'‘“ = 

w" =[/i-W^2W2i]wfi^ 



29 



APPENDIX B. EFFECTS OF CORRECTIONS 



Some selected results of applying the method are given in Figures B1 
through B4; two solutions for each of the four days, top and side views for 
each. Original track from the first member of the array pair is marked with a 
cross, and from the second member with a small circle. The (timing) 
corrected tracks are marked with connected solid lines for the former and 
dashed lines for the latter. 

The first set in Figure B1 shows no noticeable corrections, S is small. The 
second set in that figure has 5 = 3.54 and is a real data version of the situation 
depicted by Figure 3. The first set of Figure B2 appears as if one track is 
corrected more than the other. But this can be explained because the 
"stretching angles" are different. The second set of this figure shows a rather 
common condition in that the corrected track is at a shallower depth than the 
original. 

The first set in Figure B3 has something of a "showcase" flavor; the top 
view has corrected track with desirable coherence. (But a svstematic 
separation in the vertical remains.) The second set also shows tantalizing 
improvement. The first set in Figure B4 has S = -0.8, and provides an 
example of modest "contracting" of points. The second set is an extreme case 
of curved track. The top view has an excellent correction displaved, but the 
side shows that the vertical is still separated. 
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Figure Bl. Before and After Offset Correction - 3/23/89 



31 



I)cpth Cross Range Depth Cross Range 



Top View 




T=43 



Down Range xl(M 

5=3.04 Pair(UD 



-360 

-370- 

-380- 

I 

I 

-390 

2.24 



Side View 




J500, 



Top View 



5000 ' 



>500 — 



-360 












5.05 5.06 5.07 5.08 5.09 5.1 5.1 1 5.12 5.13 5.14 5.15 

Down Range xlOt 

T=50 5=4.91 Pair (5,56) 

Side View 



-370. 

-380. 

-390 







5.05 5.06 5.07 5.08 5.09 5.1 5.1 1 5.12 5.13 5.14 5.15 

Down Range xlO^ 



Figure B2. Before and After Offset Correction - 5/10/89 
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Figure B3. Before and After Offset Correction - 6/6/89 



33 



I)cpth Cross Range Depth Cross Range 



1340 



Top View 




-390 



Down Range 

T=50 5=-0.80 Pair (4,54) 
Side View 



xlO 



-400- 

-410- 

-420- 







3.86 3.88 3.9 



3.92 3.94 3.96 3.98 4 4.02 

Dowti Range xlCM 



5500 



6000- 

6500- 

7000 — 
5.705 





















5.706 



T=50 



-340 



Top View 




Down Range xlCM 

5=2.24 Pair (15,16) 

Side View 



-350 - 



-360 

-370- 

5.705 







5.706 




5.706 5.707 5.707 

Down Range 






5.708 5.709 

xl04 



Figure B4. Before and After Offset Correction - 7/21/89 
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APPENDIX C COVARIANCE MATRICES 



The covariances of the residual processes, eq. (38), are recorded here along 
with M, eq (39), which estimates M, eq. (16). The cross covariances D^y are 
converted to correlations, Rxy- by 



^xy(h/) - ^Ary(h/)/ < h/ — 1/2,3. 



The order of presentation is that of Table 1. The number of points used, T, is 
occasionally smaller than its Table 1 counterpart. This is due to outlier 
rejection based upon a visual view of the track. In three instances these 
values are zero and then, the covariances are not computed. 

In a number of instances the covariances are marked with a "c" or "cc" in 
the identifying column. This means the track was curved (c) or excessively 
cu~ved (cc) so that the straight line filter could not reasonably be assumed to 
provide a valid representation of the path. It is a curiosity that, in many of 
these instances, the cross correlations are strong and their use produces, via 
eq. (16), compensatory values. That is, the M matrices do not appear to be 
especially large. Some large cross correlations also appear in cases that pass 
the visual test for straight line track. 
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C.l, Covariances and Cross Correlations of Residuals 
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APPENDIX D. AUTOCORRELATIONS 



Examples of autocorrelation functions of residuals R(h)/R(0), (see eq. (41)) 
are presented in Figure D1 through D6. The first two are for the displacement 
process, eq. (40). The latter four are for the three components of displacement. 
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Figure Dl. Autocorrelation of noise displacement 
3/23/89 5 = 0.15 
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Figure D3. Autocorrelation of noise components 
5/10/89 6 = 1.24 
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Figure D4 Autocorrelation of noise components 
5/10/89 5 = 2.85 
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Figure D5. Autocorrelation of noise components 
6/06/89 5 = 2.48 
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APPENDIX E. COMPUTER SOURCE CODES 



The major programs needed to produce Table 1 are documented here. 
Basically there are three steps; 

(i) Extract pertinent crossover data from T-files 

(ii) Produce the covariance matrices of residuals; Dx, Dy, Dxy and M. See 
equations (38) and (39) 

(iii) Develop the estimates, eq. (9) and (10). Develop supporting 
statistics, eq (17), (18), (19) 

Also the user needs access to Table E.l, the coordinates of the various 
arrays at the Nanoose range. 

Step (i) is managed by the program KEYGATE, written in FORTRAN 77, 
Miscosoft optimizing compiler 4.01. This program reads the NUWES T-files 
and performs a series of gating operations in order to produce crossover data 
in the required seven column format. It will prompt the user for 

(a) name of the range (e.g., Nanoose) 

(b) Number of records to ignore (e.g., 31 for bypassing the DVT 
information) 

(c) The target vehicle mode (e.g. 7) 

First, it will strip out all data other than mode 2 or mode 7. Second, it 
identifies all point counts for which a position vector is available for two or 
more arrays. Third, it reads an array location file and removes data that 
cannot be reasonably identified with an overlap region. Next, it organizes the 
data by pairs of arrays. Finally, it prompts the user for his array pair and 
output filename; it selects and records all the data of the desired type. 
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Step (ii) is contained in the program PRINCOM3.M. This is a 
PCMATLAB program that performs the eigen analysis, eq. (32); develops the 
smoothed fit, eq. (37); and the covariance matrices (38) and (39). The latter is 
placed on an output file to be used in Step (iii). It also develops the 
autocorrelation sequence, e.g. (41) of the displacement process. 

A slight modification of this program can be made to develop 
autocorrelations for the three components of residuals instead of for the 
displacement process. 

Step (iii) is performed by the program TIMECOR, a FORTRAN 77 code. 
The Microsoft FORTRAN optimizing compiler 4.01 is used. The user should 
note that requirements of the three input files: 

VELOCITY.DAT is a two-column data set containing the water layer 
boundaries and the sound speeds in 25 ft increments 

TRPDOTRX.DAT is the seven-column crossover data output of KEYGATE 
to which has been appended the locations of the two arrays as the last 
record. 

MMATRIX.M is the covariance matrix, M, produced by PRINCOM.M in 
Step (ii). 

Since the ray tracing exit layer elevation angles are not contained in 
T-files it is necessary to reconstruct them in order to compute Equation (5). 
This is accomplished by the subroutine TGEN. The methodology of TGEN is 
explained in [3] under the heading of ray fitting. 

There are two output files which, taken together, contain all of the 
information indicated in Table 1. 
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TABLE E.I. COORDINATES OF THE NANOOSE ARRAYS 



Array Number 


Date of Survey 


Downrange 


Crossrange 


Depth 


0 


6/20/85 


12188.01 


-131.52 


-1295.33 


1 


6/20/85 


19463.16 


-174.99 


-1308.76 


2 


7/12/85 


26991.39 


-109.83 


-1323.25 
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1/07/88 


34505.10 


-80.76 


-1323.32 
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10/24/88 


42005.19 


-55.17 


-1318.28 
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6/20/85 


49497.00 


-25.23 


-1315.58 
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6/20/85 


56972.28 


-21.21 


-1308.50 
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7/30/85 


64680.66 


15.33 


-1353.39 
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11/16/88 


71969.73 


-29.28 


-1300.89 
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5/07/S4 


3.00 


3.00 


1.00 


10 


3/12/84 


47100.00 


-3600.00 


-1300.00 


11 


7/18/85 


23173.89 


-6488.40 


-1312.09 


12 


6/20/85 


30731.25 


-6553.05 


-1312.90 


13 


6/20/85 


38213.61 


-6640.77 


-1323.05 


14 


6/20/85 


45647.07 


-6513.18 


-1324.78 


15 


6/19/85 


53249.43 


-6354.60 


-1316.66 


16 


9/13/85 


60859.74 


-6356.07 


-1313.42 


17 


6/16/87 


68217.93 


-6524.10 


-1313.43 


54 


2/02/88 


38029.95 


5401.98 


-1212.69 


55 


6/20/85 


45645.75 


6369.66 


-1188.12 


56 


7/30/85 


53180.13 


6417.96 


-1218.&1 


57 


7/30/85 


60745.71 


6419.40 


-1088.24 
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6/20/85 


41605.14 
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-1268.23 
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-1255.35 
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7/15/80 


22119.60 
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5/04/83 


45000.00 


1500.00 


-1350.00 
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.00 


.00 



KEYGATE 



PROGRAM KEYGATE 

Program to read in raw data from Keyport hydrophone arrays, segregate it by mode, and throw out 
unusable records. The output of this program is to be read in by the program KEYMAIN. 

Modified by Colin R. Cooper 11/15/89 



INTEGERM CRT, KBD 
CHARACTER*25 D5NAME, SITNAM 
CHARACTER"9 TEMPI, TEMP2, TEMP3, TEMP4 

PARAMETER (KBD=5, CRT =6) 



Wl^ITECCRT,*) ' Please enter the name of your input file: ' 
READ(KBD,‘(A)') D5NAME 

WRITE(CRT,*) ' Please enter the name of the range configuration file:’, 
READ(KBD,’(A)‘) SITNAM 

TEMPI = TEMPI. IMF 
TEMP2 = TEMP2.TMF 
TEMPS = TEMP3.TMF 
TEMP4 = TEMP4TMF 

CALL STRMOD(CRT,KBD,TEMPl ,D5N AME) 

CALL PAIR(CRT,KBD,TEMP1,TEMP2) 

CALL RANGE(CRT,KBD,SITNAM,TEMP2,TEMP3) 

CALL PAIR2(CRT,kBD,TEMP3,TEMP4) 

CALL REC(CRT,KBD,TEMP4) 

VVRITE(CRT,*)’ Operation complete. KEYGATE terminating... ’ 

STOP 

END 



SUBROUTINE STRMOD(CRT,KBD, TEMPI, DSNAME) 

Program to strip all modes except 2‘s and 7’s from the Kevport data. (2 indicates target ship, 7 
indicates torpedo). 



CHARACTER DO*2, DSNAME^25, TEMPI "9 
INTEGER PC, ARRAY, CRT, NHEAD 

OPEN(l,nLE=DSNAME,STATUS=‘OLD') 

WRITE(CRT,*)' How manv records of header do you want to strip off the file?', 
READ(KBD,N NHEAD 
DO 11 I = 1,NHEAD 
READfLN 
11 CONTINUE 



SI 



WRITE(CRT,')' Input mode to be kept?' 

READ(KBD,‘) NUM 

OPEN(2,FILE=TEMPl,STATUS='NEW) 

10 READ(1,100,END=30,ERR=40)PC,DO,X,Y,Z,ARRAY,MODE 
IF(IX) .NE. ' •)COTO20 
IF(MODE .NE. NUM) GOTO 20 
WRITE(2,llO)PC,X,Y,Z,ARRAY,MODE 
20 CONTINUE 
GOTO 10 

40 WRITE(CRT,-)' THERE WAS AN ERROR IN THE RLE' 
50 CONTINUE 

100 FORMAT(I5,A2,1X,F7.1,2X,F7.1,2X,F7.1,30X,I2,2X,I1) 

110 FORMAT(lX,I5,2X,3n0.1,2X,I2,2X,I2) 

CLOSE (UNIT=1) 

CLOSE(UNIT=2) 

RETURN 

END 



SUBROUTINE PA1R(CRT,KBD,TEMP1,TEMP2) 

Program to pair point counst after the data has been gated by STRMOD. Second pass. 



DIMENSION X(200), Y(200), Z(200) 

INTECERM PC(200), ARRAY(200), MODE(200), CRT, HOLD 
CHARACTER'9 TEMPI, TEMP2 

OPEN(l,RLE=TEMPl,STATUS='OLD') 

OPEN(2,FILE=TEMP2,STATUS='NEW) 

HOLD = 0 
IREC = 0 
NREC = 0 
1 = 1 



Read records by two's to compare p>oint counts. 



READ(1,’,END=40,ERR=30) PC(I), X(I), Y(I), Z(I), 

+ ARRAY(l), MODE(I) 

10 READ(1,*,END=40,ERR=30) PC(1+1), X(I+1), Y(I*1), Z(I*1), 
+ ARRAY(I-l), MODE(I-rl) 

NREC = NREC + 1 

IF(PC(I+1) .EQ. HOLD) THEN 

WT?ITE(2,100) PC(I+1), X(Irl), Y(I+1), Z(I*1), 

+ ARRAY(1+1), MODE(I+l) 

HOLD=PC(Ul) 

GO TO 20 
END IF 

IF(PC(I) .EQ. PC(I+D) THEN 

WRITE(2,100) PC(I), X(I), Y(I), Z(I), ARRAY(I), MODE(I) 
v\T:n'E(2,ioo) pc(ui), xa-i), yiud, Z(I-d, 



+ ARRAY(Ul), MODE(I+l) 

HOLD = PC(I+l) 

IREC = IREC + 1 
GO TO 20 
END IF 

1F(FC(1) .NE. PC(l-fl)) THEN 
PCa) = PC(l-Kl) 

x(i) = xa+i) 

Ya) = Y(i+i) 

Z(1) = Z(1+1) 

ARRAY(l) = ARRAY(l-^l) 

MODE(1) = MODE(I-h1) 

1 = 1 
END IF 

GOTO 10 
20 1 = U1 
GO TO 10 

30 WRITE(CRT,*)‘ THERE IS AN ERROR IN THE DATA HLE IN REC/,NRE 
WRITE(CRT,*)’ ... OPERATION TERMINATING DUE TO ERROR. * 

STOP 

40 CONTINUE 

CLOSE(UNlT=l) 

CLOSE(UNIT=2) 

100 FORMAT(lX,I5,2X,F7.T2X,2F9.1,2X,12,2X,I2) 

RETURN 

ENT) 



SUBROUTINE RANGE(CRT,KBD,SITNAMTEMP2JEMP3) 

TTiis program completes the third gating of Keyport range data. It reads array location data from a site 
specific configuration file and tests to see if the data is in the valid overlap area. 



INTEGERS ARRAY, CARRAY, CRT, P 
REALM CONTIG(200,4), LX, LY, LZ, MAXTAL 
CHARACTER SITNAM"25, TEMP2"9, TEMP3*9 



Open input and output files: 



OPENU ,nLE=TEMP2,STATUS=*OLD’) 

OPEN(2,FILE=SITNAM,5TATUS=OLD’) 

OPEN(3,nLE=TEMP3,STATUS=‘NEW) 



Read site configuration into CONTIG array; 



NREC = 0 
1 = 1 

10 READ(2,’,ENT)=40,ERR=30) CONnG(U), CONnG(I,2), CONHGILO), 
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-K CONFIG(U) 

NREC = NREC -K 1 
1 = 1^1 
GOTO 10 

30 WRITE(CRT/)’ There was an error reading the config file in record^iVRB 
40 CONTINUE 
CLOSE(UNIT=2) 

NDREC = 0 



Read X, Y, Z, and ARRAY from input data file: 



45 READ(1/.ERR=70,END=80) PC, X, Y, Z, ARRAY, MODE 
DO50I = 1,NRE 
LX = COrsTFIGd,!) 

LY = CONFIG(I,2) 

LZ = CONFIG(I,3) 

CARRAY = INT(CONnG(I,4)) 

MAXVAL = 4700. 



Match array number in data file with that in config. file. If they are equal compute slant range 
distance (SkDIST): 



1F(ARRAY .EQ. CARRAY) THEN 

SRDIST = SQRT( (X - LX)^2 -k (Y - LY)**2 ^ (Z - LZ)**2 ) 
IF(MAXVAL .GE. SRDIST) THEN 

WRJTE(3,100) PC, X, Y, Z, ARRAY, MODE, SRDIST 
NDREC = NDREC + 1 
ENT) IF 
ENT IF 

50 CONTINUE 
GOTO 45 

70 WRITE(CRT,*)’ There is an error in the data file. ' 

80 CONTINUE 

CLOSE(UNTT=l) 

CLOSE(UNIT=3) 

100 FORMAT(3X,15,3(3X,F10.1),2(3X,I2)3X,F8.2) 

RETURN 

END 



SUBROUTINE PAIR2(CRT,KBD,TEMP3,TEMP4) 

Program to pair point counts after the data has been tested by RANGE. Fourth pass. 



DIMENSION X(200), Y(200), Z(200), SRDIST(200) 

INTEGERS PC(200), ARRAY(200), MODE(200), CRT, HOLD 
CHARACTER*9 TEMP3, TEMP4 
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OPEN(l,nLE=TEMP3,STATUS=‘OLD’) 

OPEN(2,nLE=TEMP4,STATUS=’NEVV) 

HOLD=0 
IREC = 0 
NREC = 0 
1 = 1 



Read records by two’s to compare point counts. 



REAEX1/,ENT)=40,ERR=30) PC(I), XdX Y(I), Z(I), 

Hh ARRAY(I), MODE(I), SRDIST(I) 

10 REAEX1,\END=40,ERR=30) PC(M), X(I+1X Yd^l), Zd+l), 

-K ARRAYd+1), MODE(Kl), SRDIST(Kl) 

NREC = NREC -f 1 

IF(PCd-l) .EQ. HOLD) THEN 

VMUTE(2,100) PC(I^l), Xd-1), Y(I^l), Z{W\), 

^ ARRAY(I+1), MODEd-1), SRDIST(I^l) 

HOLD = PC(I+l) 

GOTO 20 
END IF 

1F(PC(I).EQ. PC(kl))THEN 

V\WE(2d00) PCd), X(I), Yd), Zd), ARRAYd), MODE(I), 

^ RDlSTd) 

V\T^ITE(2,100) PCd^l), Xd^l), Y(H1), Z(I^l), 

+ ARRAY(I + 1), MODE(M), SRDIST(I+1) 

HOLD = PC(M) 

IREC = IREC ^ 1 
GOTO20 
ENT) IF 

IF(PCd) .NT. PC(K1))THEN 
PC(1) = PC(I-1) 

X(1) = X(1-1) 

Y(l) = Y(K1) 

Z(1) = Z(M) 

ARRAY(l) = ARRAY(Kl) 

MODE(l) = MODE(I-l) 

SRDlST(l) = SRDIST(Ul) 

1 = 1 

ENT) IF 

GOTO 10 
20 1 = U 1 
GOTO 10 

30 VVRITE(CRT/)' THERE 15 AN ERROR IN THE DATA RLE IN REC.',NRE 
VVRITE(CRT/)‘ ... OPEFL\TION TER.MINATING DUE TO ERROR. ’ 

STOP 

40 CONTINUE 

CL05E(UMT=1) 

CL05E(UN1T=2) 



IOC FORMATnX,I5,2X,F7.1,2X,2F9.1,2X,I2,2XJ2,2X,F8.2) 



RETURN 

END 



SUBROUTINE REC(CRT.KBDJEMP4) 

Program to produce the final gating of Keyport hydrophone array test data. 



INTEGER PC, ARRAY, PCI, PC2, Al, A2, ARRAY!, ARRAY2, CRT 
INTEGER PCAOO), ARRAYAOO), HOLD 
DIMENSION XA(IO), YA(!0), ZA(!0), SRDISTA(lO) 

CHARACTER OUTFIL*25, ANSM, TEMP4*9, TEMP!*9,RANGE*7 

RANGE=’NANOOSE‘ 

99 WRI7T(CRT,*)‘ What is the name you wish to give to ‘, 

’ your output file? 

READ(CRT,’(A)‘) OUTFIL 

WRITE(CRT,*)‘ Enter the numbers of the arrays to be paired: ' 
READ(KBD,‘^) A!, A2 

OPENd ,nLE=TEMP4,STATUS=’OLD') 
OPEN(2,FILE=TEMPl.DAT,STATUS=’OLD*) 

5 READ(1,*,END=8,ERR=8) PC, X, Y, Z, ARRAY, MODE, SRDIST 
IF((ARRAY .EQ. A!).OR.( ARRAY .EQ. A2)) THEN 
V\T^ITE(2,!00) PC, X, Y, Z, ARRAY,SRDIST 
END IF 
GOTO 5 

8 CONTINUE 
CLOSE (UNIT=1) 



Routine to pair data again: 



REWIND 2 

OPEN(4,FILE=TEMP2.DAT’,STaTUS=’OLD’) 



I = ! 

IFLAG = 0 
FIRST = 1 
HOLD = 0 
M = 0 



9 READ(2,100,END=40) PCA(I), XA(I), YA(I), ZA(I), 
+ ARRAYA(I), SRDISTA(I) 

IRHRST .EQ. DTHEN 
nRST = 0 
HOLD = PCA(I) 

ENT) IF 

IF(PCA(I) .EQ. HOLD) THEN 
I = U 1 
M = M - 1 
GOTO 9 



ELSE 



In cases where there are three or more reports for a given p>oint count, segregate by companng 
SRDIST. if there are more than 3 reports, discard all. 



IF (M .EQ. 3) THEN 

IF(ARRAYA(M) .EQ. ARRAYA(M-D) THEN 
IF (ARRAYA(M) .EQ. ARRAYA(M-2)) THEN 
GOTO 50 
ELSE 

IF (ABS(SRDISTA(M-2)-SRDISTA(M)) .LT. 

+ ABS(SRDlSTA(M-2)-SRDISTA(M-l))) THEN 

V\'RITE(4,100) PCA(M), XA(M), YA(M), ZA(M), 

+ ARRAYA(M), SRDISTA(M) 

WRITE(4,100) PCA(M-2), XA(M-2), YA(M-2), ZA(M-2), 
+ ARRAYA(M-2), SRDISTA(M-2) 

IFLAG = 1 
ELSE 

WR1TE(4,100) PCA(M-l), XA(M-l), YA(M-l), ZA(M-l), 
ARRAYA(M-l), SRDISTA(M-l) 

^^1^(4,100) PCA(M-2), XA(M-2), YA(M-2), ZA(M-2), 
+ ARRAYA(M-2), SRDISTA(M-2) 

IFLAG = 1 
END IF 
END IF 
ELSE 

IF (ARRAYA(M) .EQ. ARRAYA(M-2)) THEN 
IF (ABS(SRDISTA(M-l)-SRDISTA(M)) .LT. 

+ ABS(SRDISTA(M-l)-SRDlSTA(M-2))) THEN 

\VRITE(4,100) PCA(M), XA(M), YA(M;, ZA(M), 

+ ARRAYA(M), SRDISTA(M) 

VvT?ITE(4,100) PCA(M-l), XA(M-l), YA(M-l), ZA(M-l), 
-r ARRAYA(M-l), SRDISTA(M-l) 

IFLAG = 1 
ELSE 

V\T^ITE(4,100) PCA(M-l), XA(M-I), YA(M-I), ZA(M-l), 
^ ARRAYA(M-l), SRDISTA(M-l) 

V\'RITE(4,100) PCA(M-2), XA(M-2), YA(M-2), ZA(M-2), 
ARRAYA(M-2), SRDISTA(M-2) 

IFLAG = 1 
ENT) IF 
ELSE 

IF (ABS(SRDISTA(M)-SRDISTA(M-D) .LT. 

ABS(SRDISTA(.M)-SRDISTA(M-2))) THEN 
UT^ITE(4,100) PCA(M), XA(M), YA(M), ZA(M), 

- ARRAYA(M), SRDISTA(M) 

WRITE(4,100) PCA(M-l), XA(M-l), YA(M-l), ZA(M-l), 

ARRAYA(M-l), SRDISTA(M-l) 

IITAG = 1 
ELSE 

UT^ITE(4,100) PCA(M), XA(M), YA(M), ZA(M), 

- ARRAYA(M), SRDISTA(M) 

V\T^ITE(4,100) PCA(M-2), XA(M-2), YA(M-2), ZA(M-2), 

ARRAYA(M-2), SRDISTA(M-2) 

IFLAG = 1 
E.NT) IF 



^7 



END IF 



ELSE 

IF ((M .EQ. 2/ .AND. (ARRAYA(M) .NE. ARRAYA(M-l))) THEN 
WRITE(4,100) PCA(M), XA(M), YA(M), ZA(M),ARRAYA(M), 
+ SRDISTA(M) 

WRITE(4,100) PCA(M-l), XA(M-l), YA(M-l), ZA(M-l), 
ARRAYA(M-l), SRDISTA(M-l) 

END IF 
END IF 

50 CONTINUE 

25 CONTINUE 

IFLAG = 0 
PCA(1)=:PCA{I) 

XA(1) = XA(I) 

YA(1) = YA(I) 

ZA(1) = ZA(I) 

ARRAYA(l) = ARRAYA(I) 

SRDISTA(l) = SRDISTA(I) 

1 = 2 
M = 1 

HOLD= PCA(l) 

EX>45L = 2,10 
PCA(L) = 0 
XA(L) = 0 
YA(L) = 0 
ZA(L) = 0 
ARRAYA(L) = 0 
SRDISTA(L) = 0 
45 CONTINUE 

ENT) IF 

GOTO 9 

40 CONTINUE 
NREC = 0 
CLOSE (UN1T=4) 

OPEN(7,nLE='TEMP2.DAT,STATUS='OLD') 

0PEN(9TILE=0UTFIL,STATUS='NEVV) 

WRITE(9,300) RANGE, Al, A2 



Read in array data in two record pairs: 



10 READ(7,100,END=60,ERR=70) PCI, Xl, Yl, Zl, ARRAY!, SRDISTl 
READ(7,100,END=60,ERR=70) PC2, X2, Y2, Z2, ARRAY2, SRDIST2 

IFtARRAYl EQ. Al .AND. ARR.AY2 EQ, A2> THEN 



If arrays are in specified order (e.g. 4^): 

WRITE(9,200) PCI, XI, Yl, Zl, X2, Y2, Z2 
END IF 

1F(ARRAY1 .EQ. A2 .AND. ARRAY2 .EQ. Al) THEN 
If arrays are in reverse order (eg. 5,4): 



WTirrE(9,200) PCI, X2, Y2, Z2, XI, Yl, Zl 
END IF 



Increment record counter; 



NREC = NT^C + 1 
GOTO 10 



70 WRITECCRT/)’ There is a bad record in the file.’ 
60 CONTINUE 

CL05E(UNIT=2) 

CL05E(UN1T=7) 

CLOSE(UNHT=9) 

100 FORMAT(15,3(2X,F8.1),2X,I2,lX,F8.2) 

200 FORMAT(2X,15,lX,6(2X,Fll.D) 

300 FORK*AT(16X,A10,2X,I2,3X,12) 



WRITECCRT,*) ' Do you want to tr\- another array pair? (Y/N)’ 

READ(KBD,’(A)’) ANS 

1F( AN5 .EQ 'V .OR. ANS .EQ. >0 GO TO 99 

RETURN 

END 



PRINCOM3 

function! Au toco rx,Autocor\i=princom3(fname) 

3/13/90 

PRINCOM3.M will evaluate the principal components of the passed data file. Computes the 
autocorrelation of the displacement process. 



iname = [’d:\mdata\’ fnamc ‘.out’]; 

evaUrioad ' iname ]), 
evaUpKeyout = ' fname 
evaKl c'jcar ‘ fname]); 



59 



X=kcyout(:3:7); 

Y=keyout(:,2:4); 

(lcn,nl=sizc(keyout); 



Compute averages and covariance matices for the tracks produced by each of the arrays. 

Xav=ones(len,l)*sum(X)/len ; 

Yav=onesOen,l)*sum(Y)/len ; 



XX=X-Xav ; 
YY=Y-Yav ; 



ipc=keyout(:,l); 



CX=cov(XX); 

CY=cov(YY); 



Develop the eigen analysis. Choose largest eigen value. 

|WX,DXl=cig(CX) 

|VN%DY]=eig(CY) 

ifDX(l,l)>DX(2,2) 

if DX(1,1)>DX(3,3) 
lx=l ; 
end 

elseif DX(2,2)>DX(3,3) 
lx=2 ; 

end 

if DX(3,3)>DX(1,1) 

if DX(2,3)>DX(2,2) 
lx=3; 
end 

end 

ifDY(U)>DY(2,2) 

if DY(1,1)>DY(3,3) 
ly=l ; 
end 

elseif DY(2,2)>DY(3,3) 
ly=2 ; 

end 

ifDY(3,3)>DY(l,l) 

if DY(3;3)>DY(2,2) 
ly=3; 
end 

end 



PX=XX*WX; 

PY=\T*WY; 

ux=PX(:,lx); 

uy=PY(:,ly); 



Modify data projection onto the first principal component to account for constant speed. 



ipca\’=sum(ipc)/length(ipc); 
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mx=(sum(ux.*(ipc-ipcav)))/(sum((ipc-ipcav).*(ipc-ipcav))) ; 
my=(sum(uy.*(ipc-ipcav)))/(sum((ipc-ipcav).*(ipc-ipcav))) ; 

ax=(sum(ux)/length(ux)) - (mx^ipcav); 
ay=(sum(uy)/length(uy)) - (my'^ipcav); 

uhx=(mx*ipc)-Hax; 

uhy=(my*ipc)+ay; 

PPX=zeros(len,3); 

PPY=zeros(len,3); 

PPX(:,lx)=uhx; 

PPY(:,ly)=uhy; 



Compute the straight line tracks in the original coordinate system. 



XSL=(PPX*VVX’)+Xav; 

YSL=(PPyWY')-KYav; 



Develop the displacement process of residuals. 



distx = sqrt((X(:,l) - XSL(:,1)).'^2 -k (X(:,2) - XSL(:,2)).'"2 + (X(:,3) - XSL(:^)).^2); 
disty = sqrt((Y(:,l) - YSL(;,1)).'^2 + (Y(;,2) - YSL(:,2)).^2 + 0'(:^) - YSU:,3)).^2); 

da\"x = sum(distx)/length(distx); 
da\y = sum(disty)/length(disty); 

\'x = distx - davoc; 

\y = disty - da\y; 



Compute the auto correlations. 



Autocorx = zeros(length(vx),l); 

Autocory = zeros(length(\y),l); 

for k = 1 ;length(\y) 

for i = l:length(\y) - k -r 1 

Autocorx(k) = Autocorx(k) ^ vx(i)*\oc(i + k-1); 

Autocory(k) = AutocoryCk) + \y(i)*vy(i + k-1); 
end 

Autocorx(k) = Autocorx(k)/sum(\^.^2); 

AutocoryCk) = Autocory(k)/sum(\y.^2); 

end 

clg 

subplot(2ll) 

plot(Y(;,l),Y(;,2);or',YSL(:,l),YSL(:,2);-r*,X(;,]),X(:,2)/xg‘,XSL(:,l),XSL(:,2);-g’) 
titleCifname ‘ - Corrected Tracks with Principal Components')) 
xlabcK'Dow’n Range')ylabelfCross Range’) 
subplot(212) 

plot(YC;,l),Y(:,3);or',YSL(:,l),YSL(:,3);-r,X(;,l),XC:,3j;xg\ X5L(:,l),XSL(;,3);-g') 
titlerCorrected Tracks with Pnncipal Components’) 
xlabelTDown Rangc‘),ylabel(’Depth‘) 
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pause 

clS 

h = 0:length(Autocorx)-l; 

subplot(211),plot(h,Autocory/or ,h,AutocoTy/-r') 
titledfnamc ' - Autocorrelation of Distance, Y Array'D.grid 
subplot(212),plot(h,Autocorx/or',h,Autocorx/-r’) 
titkK’Autocorrelation of Distance, X Array’), grid 

plot(h(l :27),Autocorx(l :27),‘or',h(l :27),Autocorx(l ;27),’-r') 
titleCAutocorrelation of Residual Distancc'),grid 
xlabeK’Point Count Lag’) 



TIMCOR 



PROGRAM TIMECOR 



06/06/90 

This program estimates the timing synchronization offset and drift parameters based upon cross-over 
data from T-files. Since exit angles and transit times are not recorded to these files they must be re- 
estimated. Location information for each of the involved arrays is also required. 



Tliis program was compiled using; 

Microsoft FORTFLAN OPTIMIZING COMPILER ver. 4.01 



This files must be able to access the following files 



VELOCITY.DAT 

TRPDOTRX.DAT 

MMATRIX.M 

OUTPUTl.COV 

OLTPUT2.COV 



- Sound velocity vs. depth data 

- Torpedo tracking data. 

- Correlation matrices for data. 

- Output file containing delta and tnot 

- Output file cuntdining cross caireldtion and M matrices. 



Users Notes 

- VELOCITY.DAT file contains layer boundaries and sound velocities for 25 ft. depth increments 

- TRPDOTRX.DAT file contains point counts and three position components for each of the two 

contributing arrays. The position data from the low'er numbered array is columns two through 
four. The last record contains the coordinate position of the two arrays.. 

- MMATRIX.M file contains the set of covariance matrices produced by AUTCXZOR3.M. 

- The two output files contain all the information appearing in Table 1. 



DIMENSION ircnso) 

CHARACTER*30 LINE 

REAL-^ Y1(150,2),Y2(150,2), Y3(1S0,2),YC1 (150,2), YC2(1S0,2) 
REALMS LL(55),G(55),W(55),A2,P1,P2,V0,V1,TT1ME(150,2) 
REALMS THETAn50,2),DEPTH(55),PHl(150,2),Vn50,2),DZ 
RE AL*8 B1 (1 50,2),B2( 1 50,2), B3( 1 50,2), YC3(1 50,2) 

REAL-^S DEN,GG,SSB,SSR,TD,TM,SIG,DD,TEMP 
REAL’^S ME(2,GEQ(150),DELTA,T1,TNOT,D1,D2,MSR,MSB 
REALMS DUM1,DUM2,DUM3,DUM4,DLJM5,DUM6,DUM7 
REALMS M(3,3),C(150,3),VDEL,SDEL,VM,SDVM,COVDM 
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REAL’S DELTAN,MEQN,R,CX(3,3),CY(3^),CXY(3^),RXY(3,3),T 
INTEGER’4 TEE 



Read in the data from the data files. 



OPEN(UNIT=2,FILE=*VELOCITY.DAT,STATUS=*OLD‘) 

OPEN(UNIT=7,nLE=TRPDOTRX.DAr,STATUS=’OLD’) 

OPEN(UMTT=lO,nLE=*MMATRIX.M‘,STATUS=’OLD') 

OPEN(UNIT=ll,nLE='OLrrPUTl.COV*,STATUS=’OLD') 

OPEN(UNIT=12,nLE=*OUTPUT2.COV,STATUS=’OLD’) 

RE AEX10,1 20)CX(T1 ),CX(2,1 ),CX(3,1 XCX(L2),CX(2,2L 

• CX(3,2LCX(1 ,3),CX(2,3),CX(3,3) 
READ(10T20)CY(T1),CY'(2,1),CY(3,1XCY(1,2),CY(2,2L 

’ CY'(3,2),CY(1 ,3XC Y(2,3),CY(3,3) 

RE AEXl 0,1 20)CXY(1 ,1 ),CXY(2,1 ),CX Y(3,l ),CX^'(1 ,2),CXY(2,2), 

• CXY(3,2),CXY(1,3),CX^'(2,3),CXY(3,3) 

REAEX10,124)T 

TEE = IDINTCT) 

IFCTEE.EQ.O) THEN 
D02I = 1,3 

RXY(1,I) = 0.0D0 
RXY(2,I) = O.OEX) 

RXY(3,I) = O.ODO 

2 CONTINUE 
ELSE 

0041 = 1,3 
EX)3J = 1,3 

RX'Y(j,l)=CXT(j,I)/DSQRT(CX(U)’CY(JJ)) 

M(],I) = CXO J) CY(j,I) . CXT(I,J) - CXT(U) 

3 CONTINUE 

4 CONTI aNUE 
ENDIF 

5 REAEK12;(A)',END=7)LINE 
GOTO 5 

7 BACKSPACES 
WRITE(12,122)TEE 
EX)8I = 1,3 

8 \VRITE(12,121>CX(L1),CX(L2),CX(1,3),CY(I,1),CY(I,2), 

• CY(I,3),RXY(I,1),RXY(I,2),RXY(I,3),M(I,1),M(I,2),M(I,3) 
WRITE(12,123) 

1=1 

10 READ(7,*)IPC(I;,Y1(U),Y2(L1),Y3(I,1), 

IF(IPC(I).EQ.999) GOTO 15 
I=H1 
GOTO 10 
15 LEN=M 

120 FORMAT(lX,9(E15.4,lX)) 
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121 FORMAT(5X,I2,5X,i •,25X,1 ',2X,1 ',25X/I ■,2X/I ',25X,1 ' 

2X/I ',25X;I •) 

122 FORMAT(12X;rMX,3(F7.2,lX);l,2X;r',lX,3(F7.2,lX);T, 

• 2X;rMX,3(F7.2,lX)/T,2X,fMX,3(F7.2,lX)/l) 

123 F0RMAT(12X,1',25X/J’,2X;L',25X;J',2X;L’,25X;J' 

' 2X,1',25X;J') 

124 FORMAT(lX,El5.4) 



Read the VELOCrTY.DAT file and prepare for isogradient raytracing. 



DZ=25 

1=1 

25 READ(2/,END=30) TEMP,VV(I) 
1 = 1+1 
GOTO 75 
30 CONTINUE 

GG=(W(M)-VV(I.7))/6. 

DO 35 J = 135 

VV(J)=VV(JT)+GG 
35 CONTINUE 
LL(1)=125 
DEPTH(1)=0 
00 401=235 

LL(I)=LL(I-1KDZ 
DEPTH(I)=DEPTH(I-1KDZ 
40 CONTINUE 
00 431=235 

G(I-1)=C\^V(I)-W(M))/DZ 
43 CONTINUE 

G(55)=G(54)+GG/DZ 
00 451=132 
45 CONTINUE 
OO 80 J=U 
N=15 
IFLAG=0 
00 70 T=1,LEN 



Set variables the the call to TGEN subroutine. TGEN returns the time and elevation. 



A2=Y3(LEN+1J) 

P1=SQRT((Y1(ITJ)-Y1(LEN+1J))*'^2+(Y2(ITJ)-Y2(LEN+1J))^^2) 

P2=-Y3(TJ) 

IF(T.NT.l) THETA(ITJ)=THETA(IT-1J) 
CALLTGEN(LUG,VV,A2,PUP2THETA(ITJ),TTIME(TJ), 
V0AT,DEPTH,IFLAG) 

IFLAG=1 



Find the azimuth angles from each array. 



PHI(IT,J)=DASIN((Y2(TJ)-Y2(LEN+1J))/SQRT((Y2(TJ)- 
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Y2(LEN-hlJ))'^*2 (Y1(ITJ)-Y1(LEN+1J))^2)) 

1F((Y1(ITJ)-Y1(LEN-h1J)).LT.O) PHI(TTJ)=3.1 41 59265359 
-PHKITJ) 



Locate the layer containing the source and set it's velocity. 



50 IF(-Y3(ITJ).LE.DEPTH(M).AMD.-Y3(ITJ).GT.DEPTH(N-1)) GOTO 60 

IF(-Y3(ITJ).GT.DEPTH(N)) THEN 
* N=NV1 

ELSE 
N=N-1 

IF(N.LE.1)THEN 

N=1 

GOTO 60 
ENDIF 
ENDIF 
GOTO 50 

60 V(ITJ)=VV(N) 



Calculate the B(T) adjustments (Sperical Coordinates : Equtaion (5)) 



Bl(ITJ)=V(ITJ)*IX:OSaHETAaTJ))*EX:C>S(PHI(ITJ)) 

B2(ITJ)=V(lTJ)*DCOS(THETA(ITJ))*DSIN(PHI(ITJ)) 

B3(ITJ)=V(ITJ)*D51N(THETA(ITJ)) 

70 CONTINUE 
80 CONTINUE 
DO90I=LLEN 
90 CONTINUE 



Calculate DELTA and TNOT. 



D2=0 

T1=0 

DEN=0 

DO 200 IT^ULEN 

D1=(((B1(TT)-B1(T,2))*(B1(ITT)-B1(IT,2)))+ 

• ((B2(ITT)-B2(T,2))*(B2(IT,1)-B2(T^)))^ 

• ((B3(TT )-B3(T,2))*(B3(IT,1 )-B3(IT,2)))) 
DEN=DEN-D1 
Tl=TK(IPC(ITrDl) 

D2=D2-(((B1 (T,l )-Bl (IT,2))*(Y1 (ITT )-Yl (IT,2))H 

• ((B2(rTT)-B2(T,2))*(Y2(ITT).Y2(IT,2)))+ 

• (B3(ITT)-B3(IT,2))*CY3(ITT)-Y3(IT,2)))) 

C(TT) = B1(IT,2)-B1(ITT) 

C(T,2) = B2(T,2)-B2(ITT) 

C(IT,3) = B3(T,2) - B3(IT,1) 

200 CONTINUE 

TD= DSQRT(DEN) 
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DEN=DEN/LEN 

TNOT=m/LEN)/DEN 

DELTA=-(D2/LEM)/DEN 



Calculate m and g. 



210 



220 

230 



D1=O.ODO 
D2=O.ODO 
DO 210 rr=i,LEN 

DD=(((Bl(rr,l)-Bl(IT,2))*(Bl(IT,l)-Bl(IT,2)))+ 
(B2(IT,1 )-B2(IT,2))‘(B2(IT,l )-B2(rr,2)))+ 
((B3(IT,l)-B3(IT,2))'(B3(IT,l)-B3(IT,2)))r 
(IPC(IT)-TNOT)’(IPC(IT)-TNOT» 

D1 = D1 + DD 

D2=D2+(((Bl(rr,l)-Bl(IT,2))*(Yl(IT,l)-Yl(rr,2)))+ 

((B2(IT,1)-B2(IT,2))*(Y2(IT,1)-Y2(IT,2)))+ 

(B3(lT,l)-B3(lT,2)r(Y3(IT,l)-Y3(IT,2))»* 

(IPC(IT)-TNOT) 

CONTINUE 

MEQ=-(D2/LEN)/(D1 /LENO 
TM = Dl 



DO230J=l,2 

DO 220 rr=l,LEN 

YCl(IT,J)=Yl(ITJ)+(Bl(rrj)*(DELTA+(MEQ’(IPC(nr)-TNOT)))) 

YC2(IT,J)=Y2(ITJ)+(B2(ITJ)*(DELTA+(MEQ*(IPC(IT)-TNOT)))) 

YC3(IT,J)=Y3(ITJ)+(B3(rrj)‘(DELTA+(MEQ‘(IPC(IT)-TNOT)))) 

CEQ(IT)=DELTA+(MEQ‘(IPCaT)-TNOT)) 

CONTINUE 

CONTINUE 



Create the output files. 



SSR = O.OEX) 

SSB = O.ODO 
IX) 231 IT = 1,LEN 

SSR = SSR + (YC1(T,1) - YC1(IT,2))"»2 + (YC2(IT,1) - 

• YC2(IT,2))"2 * (YC3(IT,1) - YC3(IT,2))‘*2 
D1=(((B1(IT,1)-B1(IT,2))‘(B1(IT,1)-B1(IT,2)))+ 

• ((B2(IT,1)-B2(IT,2))*(B2(IT,1)-B2(IT,2)))+ 

• ((B3(IT,1)-B3(IT,2))’(B3(IT,1 )-B3(IT,2)))) 

SSB = SSB + DrGEQ(IT)**2 

231 CONTINUE 

MSR = SSR/(LEN - 2) 

MSB = SSB/2 
SIC = DSQRT(MSR) 

TD = TD*DELTA/SIG 
TM = MEQ*DSQRT(TM)/SIG 



Calculate VDEL,SDEL,VM,SDVM ( variances and standard deviations ). 



66 



DUM2 = O.ODO 
DUM3 = O.ODO 
DUM5 = O.ODO 
DUM6 = O.ODO 
D0 235I=i;5 
D0 233J = U 
DUMl = O.ODO 
DUM4 = O.ODO 
D0232IT-1,LEN 

DUMl = DUMl + (CaT,I)*C(ITJ)) 

DUM4 = DUM4+((IPC(IT)-TNOT)**2 *C(IT,I)*C(ITJ)) 

232 CONTINUE 

DUM2 = DUM2 + (M(IJ) * DUMl) 

DUM5 = DUM5 -e (M(IJ) * DUM4) 

233 CONTINUE 

D0 234IT= ULEN 

DUM3 = DUM3 + (C(IT,I)**2) 

DUM6 = DUM6 + ((IPC(IT) - TNOT)**2 * C(rU)*^2) 

234 CONTINUE 

235 CONTINUE 

VDEL = DUM2/DUM3**2 
SDEL = DSQRT(VDEL) 

VM = DUM5/DUM6**2 
SDVM = DSQRT(VM) 

DUM2 = O.ODO 
EX) 238 IT = l.LEN 
DUMl = O.ODO 
D0 237UU3 
DO 236 J = 1,3 

DUMl = DUMl + (C(IT,I)*C(ITJ)‘M(IJ)) 

236 CONTINUE 

237 CONTINUE 

DUM2 = DUM2 ^ ((IPC(IT) - TNOT) * DUMl) 

238 CONTINUE 

COVDM = DUM2/(DUM3 * DUM6) 

270 REAEXli;(A)’,END=280)LINE 
GOTO 270 

280 BACKSPACE 11 

IF(SDEL.EQ.O.ODO) THEN 
DELTAN = SDEL 
MEQN = SDEL 
R = SDEL 
GOTO 285 
ENDIF 

DELTAN = DABS(DELTA)/SDEL 
MEQN = DABS(MEQ)/SDVM 
R = COVDM /(SDEL^SDVM) 

285 DELTA = DELTA * lOOO.ODO 
MEQ = MEQM000.0EX) 

SDEL = SDEL MOOO.ODO 
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SDVM = SDVM • lOOO.ODO 

WRITE(11,350)DELTA,SDEL,DELTAN,MEQ,SDVM,MEQN,R,TNOT,LEN 

350 FORMAT(F5.2,2X,F5.4,2X,F6.2,2X,F7.5,2X,F7.6,2X,2(F72,2X), 

• F8.2,2X,I2) 

WRITEC,’)' PROGRAM COMPLETED !' 

END 



SUBROUTINE TGEN(LUG,VV,A2,Pl,P2,ANGLE,TIME, 
V0,V1,DEPTH,IFLAG) 



1 

1 TGEN generates transit time and elevation angle at a target if given the horizontal range, the depth of 

Ithe sensor and the target, the layer boundanes and the gradients. Isogradiant raytracing is used. 


Calling Arguments are as follows. 


LL 


- An array containing the layer midpoints. 


G 


- An array containing the gradients for each layer. 


1 V V 


- An array containing the velocity at each layer. 


A2 


- The depth of the sensor ( positive down ). 


PI 


- Range of the target ( horizontal down ). 


P2 


- Depth of the target ( positive down ). 


VO, VI 


- The values for a straight line single layer regression of depth vs. velocity. 


DEPTH 


- An array containing the depth of each layer. 


Return arguments are as toliows: 


ANGLE 


- The final angle at the target. 


TIME 


- The time of transit. 


User Notes. 




All floating point numbers are defines . REAL'S 


All times 


are in seconds, and all angles in radians. 



DIMENSION L(55),G(55),V(55),LL(55),W(55) 

DIMENSION TH(55)T(55),VZ(55) 

DIMENSION C2(55XTT(55),DEPTH(55) 

REALMS L,G,LL,VVJH,VZ,C2,TTJ,R0,A1 
REAL*8 A2,P1,P2,C1,C22,THETA,VM 

REALMS THETAZ,RV,R,TIME,ANGLE,EP,DZ, DEPTH, V0,V1,V 



Initialization: Set the value for DZ, the layer thickness. The sensor is assumed to be at 
RANGE 0. Determine the values for J, which is 1+Number of layers less than or equal to the 
sensor depth, and 1, which is the number of layers less than or equal to the torpedo depHh. 
Redefine the endpoints of those layers locally to be the depths of the torpedo and sensor. 
Define local values for the LL and arrays. 
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EP=0.1D^ 

DZ=25.0D0 

A1=0.0EX) 

J=1 

1=0 

D010K=135 

UK)=LUK) 

V(K)=VV(K) 

IF(DEPTH(K),LE.A2)J=)-h 1 
IF(DEPTH(K).LE.P2) 1=1^1 
10 CONTINUE 
IF(I.LE.O) 1=1 
N=Uj-I 

V(l)=V(lKG(inP2-UI)) 

V(J)=V(J-lHG(J-im2-L(J-D) 

L(I)=P2 

IF(A2.GT.L(J-1))G0)=G(J-1) 

L(J)=A2 

R0=P1-A1 

V1=(V(J)-V(1))/(A2-P2) 
V0=V(J) - V1*A2 



Calculate an initial estimate for the angle and time using a single layer approximation. 



C22=-V0/V1 

Cl=((0.5D0r(Pl-AD) 

C1=CU((0.5DO)"(L(I)-LO))*(L(IHL(J)-2.0DO* 
* C22)/(P1-AD) 

IF(IFLAG.GT.O) GOTO 43 
THETA=DATAN((A1-C1)/(A2<22)) 
WRITEr/)’ DEHNE THETA AGAIN* 
GOTO 49 

48 THETA=ANGLE 

49 CONTINUE 



Use the angle THETA to raytrace back through all the layers. First, use the ray invariant (RV) 
and the velocity to calculate the entrance andgle at each layer.. 



R=A1 



Find the maximum sound speed 



VM = V(I) 

DO 495 K = I J 

IF (V(K).GT.VM) VM = V(K) 
495 CONTINUE 

50 RV=DCOSTHETA)/V(J) 
DO60K=IJ 

TH(K)=DACOS(RV*V(K)) 
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T(K)=DTANCTH(K)) 

VZ(K)=V(K)-L(K)*G(K) 

C2(K)=-VZ(K)/G(K) 

60 CONTINUE 



Using the angle just calculated, iterate backwords through the layers from sensor to target to get 
the horizontal range. Stop at the depth of the target. 



R=0.0D0 

DO70K=J,I+l,-l 

Cl=R-T(Kr(L(K)-C2(K-D) 
R=CUT(K-1)*(L(K-1)-C2(K-D) 
70 CONTINUE 



Test if the value for the range is within torerance. If not, redefine THETA, the initial angle, 
and raytrace again. If wthin tolerance, calculate the time of travel based on THETA, and return. 



EP=0.1D^ 

IF((DABS(R-P1)).LE.EP) GOTO 100 
THETAZ=THETA 

THETA=DATAN(DTAN(THETAZ)*(R-A1)/R0) 

GOTO 50 

100 TT(J)=DLOG((1.0D0+DSIN(TH(J)))/(DCOSaH0)))) 
TIME=0.0D0 
EXD110K=J-1,1,-1 

TT(K)-DLOG((1.0[X)-fDSlN(TH(K)))/(DCOS(TH(K)))) 
TIME=TIME^(TT(K)-TT(K^1 ))/G(K) 

110 CONTINUE 
ANGLE=THETA 
RETURN 
ENT) 
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